Load data

Here we load the list of species that we sampled eye size for - the subset that we want to trim the phylogeny to match.

# Load salamander data
caudata_raw <- read.csv("../Data/Tidy/salamanders_tidy.csv", header=TRUE, na.strings=c("", "NA", " ")) %>%
  # Remove the gilled version of Pleurodeles_waltl
  filter(Genus_Species != "Pleurodeles_waltl_gilled") %>%
  # Remove blind salamanders from dataset
  filter(Genus_Species != "Eurycea_rathbuni" & Genus_Species != "Proteus_anguinus")

# Check out structure
str(caudata_raw)

# Tidy adult data
tip_species <- caudata_raw %>% 
  dplyr::select(Genus_Species, Family) %>%
  unique()

# Check for duplicates of species
n_occur <- data.frame(table(tip_species$Genus_Species))

Load published phylogeny

Here we import an amphibian tree published by Jetz and Pyron (2019). It’s important to note that this tree includes branches supported by molecular data as well as branches that are grafted on based on taxonomy.

We check whether the tree is rooted (it is) and whether it is reading as ultrametric (it should be, but is not, so we force it using force.ultrametric) and then look at the full tree.

# Import ultrametric phylogeny (modified) from Jetz and Pyron 2019
tree_orig <- read.tree(file = "../Data/Raw/amph_shl_new_Consensus_7238.tre") #reads tree

# Check whether tree is rooted
is.rooted(tree_orig) 
## [1] TRUE
# Check whether tree is dichotomous (no polytomies)
is.binary(tree_orig) 
## [1] FALSE
# Check that tree is ultrametric
is.ultrametric(tree_orig)
## [1] FALSE
# Force ultrametric
tree <- force.ultrametric(tree_orig)
## ***************************************************************
## *                          Note:                              *
## *    force.ultrametric does not include a formal method to    *
## *    ultrametricize a tree & should only be used to coerce    *
## *   a phylogeny that fails is.ultrametric due to rounding --  *
## *    not as a substitute for formal rate-smoothing methods.   *
## ***************************************************************
# Check that tree is ultrametric
is.ultrametric(tree)
## [1] TRUE
# Plot tree
plot.phylo(tree, show.tip.label = F)

Match names in tree and phylogeny

First, we find species that are in our dataset that aren’t found as exact matches in the phylogeny tip labels.

# Find species in dataset that don't exist as phylogeny tip label
tip_species[which(!tip_species$Genus_Species %in% as.vector(tree$tip.label)), ]

There is 1 species that doesn’t appear in the phylogeny, Aneides niger. We can look for a synonym or close relative in the phylogeny tips.

# Look for close relative in phylogeny tips
phylo.tips <- as.data.frame(tree$tip.label)

After looking at Frost (2022) and searching through the phylogeny tips, Aenides niger is not found in the tree, but it previously was a subspecies of Aneides flavipunctatus, which is present in the tree. We also have that species in our dataset, so we will graft on a branch that is sister to A. flavipunctatus for A. niger.

# Copy tree for editing
tree_graft <- tree

# Find most recent common ancestor of Aneides flavipunctatus and sister clade of Aneides_vagrans + Aneides_ferreus
bind.node <- findMRCA(tree_graft, 
                      tips = c("Aneides_flavipunctatus", "Aneides_vagrans", "Aneides_ferreus"), 
                      type = "node")

# Add A. niger to that node
tree_graft2 <- bind.tip(tree = tree_graft, 
                        tip.label = "Aneides_niger", 
                        where = bind.node)

Now we can check that our names in the dataset are matching those in the phylogeny. Now that we’ve consolidated subspecies into species, we have 155 species with data that we want to match the tree.

# Find difference between tree and dataset
rownames(tip_species) <- tip_species$Genus_Species
overlap <- name.check(tree_graft2, tip_species)

# Species in data that aren't matching tree
overlap$data_not_tree

All the species are now matching! Now we prune the tree to match the species in our dataset.

# Make list of taxa to drop (in tree but not in dataset)
drops <- setdiff(tree_graft2$tip.label, tip_species$Genus_Species)

# Drop unwanted tips from phylogeny
tree.pruned <- drop.tip(phy = tree_graft2, tip = drops) 

# See which tips you've kept in your phylogeny
tree.pruned$tip.label
  
# Check that phylogeny tips and data match exactly (if they match will return "OK")
name.check(tree.pruned, names)

The pruned tree tip labels and the dataset species names match exactly. We can plot the tree to check that it looks good.

# Plot pruned tree with polytomies
plot.phylo(tree.pruned, show.tip.label = TRUE)

We also check that the tree is still rooted and ultrametric.

# Confirm that tree is rooted
is.rooted(tree.pruned) 

# Check that tree is still ultrametric
is.ultrametric(tree.pruned)

Finally, because we added a polytomy our tree will not be binary. For analyses, we will need it to be, so we can randomly resolve polytomies by adding in tiny differences so that we can use statistical methods like PGLS.

# Test whether tree is dichotomous (shouldn't be yet)
is.binary(tree.pruned)

# Randomly resolve polytomies to make tree dichotomous using multi2di function in ape
dich.tree <- multi2di(tree.pruned)

# Check that tree is now binary
is.binary(dich.tree)

The tree is now dichotomous. We will take one final look at the tree to check that it looks ok and then export as a nexus file for use in our analyses.

# Plot final tree
plot.phylo(dich.tree, show.tip.label = TRUE)

Looks good! Now we export.

# Export final tree
write.nexus(dich.tree, file = "../Data/Tidy/caudata-tree.nex")